float4 product_q(const float4 sv, const float4 tw);
float4 rotate_q(const float4 q, const float4 p);
float nextIntersection(float4 origin, float4 dir, float4 A, float4 B);


kernel void castRay(
		global int* nodePool,
		global int* dataPool,
		read_only image3d_t brickPool,
		write_only image2d_t image,
		constant float* cameraParams) {

// Ray Casting algorithm
// *Set up variables
	constant sampler_t brickSampler;

	local float4 cameraPosition, cameraOrientation, localFwd, localLeft, localDown, nearPlaneCenter;
	local float farPlaneDistance, nearPlaneDistance, fieldOfView, nearPlaneSize;
	local int imageWidth, imageHeight;

	private register int lid = get_local_id(0);
	private register int gid = get_global_id(0);

	if (gid == 0) {
		brickSampler = CLK_FILTER_NEAREST | CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE;
	}

	if (lid == 0) {
		// cameraParams:
		// 0-3: cameraPosition
		// 4-7: cameraOrientation
		// 8: nearPlaneDistance
		// 9: farPlaneDistance
		//10: fieldOfView
		//11: imageWidth
		//12: imageHeight
		cameraPosition = vload4(0, cameraParams);
		cameraOrientation = vload4(4, cameraParams);

		float4 local;
		local.x = 0;
		local.y = 0;
		local.z = 1;
		local.w = 0;
		localFwd = rotate_q(cameraOrientation, local);

		local.z = 0;
		local.x = -1;
		localLeft = rotate_q(cameraOrientation, local);

		local.x = 0;
		local.y = -1;
		localDown = rotate_q(cameraOrientation, local);

		nearPlaneDistance = cameraParams[8];
		farPlaneDistance = cameraParams[9];
		fieldOfView = cameraParams[10];

		imageWidth = (int)(cameraParams[11]);
		imageHeight = (int)(cameraParams[12]);

		nearPlaneCenter = cameraPosition + (localFwd * nearPlaneDistance);
		nearPlaneSize = 2.0f * nearPlaneDistance * native_tan(fieldOfView / 2.0f);
	}
	barrier(CLK_LOCAL_MEM_FENCE);

// *Initialize ray at near plane
//		-determine pixel location
	//		-compute ray direction
	private int2 indices;
	indices.x = id % width;
	indices.y = (int)((float)id / (float)width);

	private float4 pixelLocation = nearPlaneCenter;

	register float ri, rj; // i and j relative to the center of the near plane
	ri = indices.x - (imageWidth * 0.5f);
	rj = indices.y - (imageHeight * 0.5f);

	// pixel's distance from the center
	ri = native_divide(ri, imageWidth) * nearPlaneSize * 0.5f;
	rj = native_divide(rj, imageHeight) * ((nearPlaneSize * native_divide(imageHeight, imageWidth)) * 0.5f);

	pixelLocation -= (ri * localLeft);
	pixelLocation -= (rj * localDown);

	private float4 direction = (pixelLocation - cameraPosition);
	direction = fast_normalize(direction);

// *Compute first sample location at either scene volume or far plane intersection

	register float isect = nextIntersection(pixelLocation, direction, )

// *Evaluate sample

//		-if we're at the far plane, we're done

//		-otherwise, descend through the octree, beginning with the root node

//			-does the current node have children?

//				-if yes:

//					-compute the needed child node index

//					-set as current node

//					-descend from there

//				-otherwise

//				-is this node a brick?

//					-if no, adjust and accumulate its contribution

//					-if yes, compute brick texture coordinates

//						-accumulate contribution from the brick texture

//		-compute next sample location, repeat

// *Store the final color value
}

float4 product_q(const float4 sv, const float4 tw) {
	// (s + v)(t + w) = (st - v.w) + (sw + tv + v x w)

	float4 v, w;
	v.x = sv.x;
	v.y = sv.y;
	v.z = sv.z;
	v.w = 0;
	w.x = tw.x;
	w.y = tw.y;
	w.z = tw.z;
	w.w = 0;
	float qs = ((sv.w * tw.w) - dot(v, w));
	float4 sw, tv;
	sw = tw;
	sw.w = sv.w;
	tv = sv;
	tv.w = tw.w;
	float4 qv = (sw + tv + cross(v, w));
	qv.w = qs;

	return qv;
}

float4 rotate_q(const float4 q, const float4 p) {
	float4 qs;
	qs.x = 0 - q.x;
	qs.y = 0 - q.y;
	qs.z = 0 - q.z;
	qs.w = q.w;

	float4 pp = product_q(product_q(q, p), qs);
	pp.w = 0;

	return pp;
}

float nextIntersection(float4 origin, float4 dir, float4 A, float4 B) {
	float4 invdir = native_recip(dir);
	float isect = -1.0f;

	float tx1 = (A.x - origin.x) * invdir.x;
	float tx2 = (B.x - origin.x) * invdir.x;

	float tmin = min(tx1, tx2);
	float tmax = max(tx1, tx2);

	float ty1 = (A.y - origin.y) * invdir.y;
	float ty2 = (B.y - origin.y) * invdir.y;

	tmin = max(tmin, min(ty1, ty2));
	tmax = min(tmax, max(ty1, ty2));

	if (tmax >= tmin) {

		float tz1 = (A.z - origin.z) * invdir.z;
		float tz2 = (B.z - origin.z) * invdir.z;

		tmin = max(tmin, min(tz1, tz2));
		tmax = min(tmax, max(tz1, tz2));

		isect = (tmin > 0.0f) ? tmin : tmax;
	}

	return isect;
}
